# cython: profile=True

"""
solves underdetermined system of equations with minumum L_infinity (Chebyshev) norm

Reference:
Abdelmalek, N.N., Minimum L_infinity solution of underdetermined
systems of linear equations, Journal of Approximation Theory,
20(1977)57-69
"""

from numpy import *

def test():
	"""
	Aeq = zeros((2,4))
	Aeq[0] = [7,-4,5,3]
	Aeq[1] = [-2,1,5,4]
	beq = array([-30,15])
	"""
	beq = array([2,2,1,-4])
	Aeq = zeros((4,6))
	Aeq[0] = [2 ,1,4, 0,-3 ,1]
	Aeq[1] = [2 ,-1,4, 0,-3 ,1]
	Aeq[2] = [5 , 1, -3, 1,-2 ,0]
	Aeq[3] = [1 ,-2,-1,5, 1, -4]
	n,m = Aeq.shape
	a = solver(Aeq,beq)
	print 'solution is', a
	print mat(Aeq)*(mat(a[:m]).transpose())

	
def solver(Aeq,beq):
	""" main solver function """
	smallno = 1e-8
	n,m = Aeq.shape
	eqc = zeros(2*n)  # column nos corresponding to given equations in dual problem
	eqc[:n] = linspace(1,n,n)
	eqc[n:] = linspace(m+n+1,m+n+n,n)	
	nullvec = zeros(m+1)
	table = zeros((m+1,n+1)) # reduced simplex tableau, last column contains basic feasible solution
	
	# Doing initialization
	table[:m,:n] = Aeq.transpose()
	table[m,n] = 1
	g = zeros(2*(m+n))  # contains cost coefficients of dual problem
	g[:n] = beq
	g[m+n:m+2*n] = -1*beq
	td = linspace(1,n,n) # contains position of column of condensed simplex tableau in dual problem matrix
	bd = ones(m+1)  # contains columns of basis matrix
	
	for i in range(m):
		if (table[i,0]>0):
			table[i,:n] = table[i,:n]*-1
			bd[i] = i+n+(m+n)+1
		else:			
			bd[i] = i+n+1			 
	bd[-1] = 1	
	for i in range(n):
		table[m,i] = -sum(table[:,i])
	
	gj(table,m,0,m,n) # gauss-jordan elimination	

	# calculating the objective function
	gb = zeros(m+1)  # cost coefficient of columns in basis matrix
	for i in range(m+1):
		gb[i] = g[bd[i]-1]
	
	z = dot(gb,table[:,-1]) # objective function, L_infi norm in this case

	# making objective function positive	
	if (z<0):
		for i in range(n):
			td[i] = td[i]+m+n

		for i in range(m+1):
		
			if (bd[i]<=m+n):
				bd[i] = bd[i]+m+n
			else:
				bd[i] = bd[i]-m-n			

	mcost,track,z = margcost(table,m,n,td,bd,gb,g,eqc)
	
	ijk = 0
	while(min(mcost)<0):	
		# calculating the column which inters the basis matrix
		argminimum = argmin(mcost)
		
		if (argminimum > n-2):
			column = track[1,argminimum-(n-1)]
			tdi = track[0,argminimum-(n-1)]
			if(tdi in eqc):
				table[:,column] = -1*table[:,column]
			else:
				table[:,column] = 2*table[:,-1] - table[:,column]
			if (tdi<=m+n):
				td[column] = tdi+m+n
			else:
				td[column] = tdi -m-n
			
		else:
			column = track[1,argminimum]

		row = -100
		bigno = 1e+20
				
		for i in range (m+1):
			if (table[i,column]>smallno):
				tmp = table[i,-1]/table[i,column]
				if (tmp<=bigno):
					bigno = tmp
					row = i

		if (row ==-100 or ijk > 25):  # brakeing condition
			return NaN*ones(m+1)	
		for i in range(n):
			if (td[i] in bd):
				td[i] = bd[row]
				table[:,i] = nullvec
				table[row,i] = 1
		bd[row] = td[column]
		
		gj(table,row,column,m,n)
		bd[row] = td[column]
		mcost,track,z = margcost(table,m,n,td,bd,gb,g,eqc)
		ijk = ijk+1
		
	a = getresult(bd,track,mcost,z,m,n) # solution vector
	return a


def margcost(table,m,n,td,bd,gb,g,eqc):
	"""calculates marginal cost"""
	
	for i in range(m+1):
		gb[i] = g[bd[i]-1]
	z = dot(gb,table[:,-1])
	mcost = zeros(2*(n-1)) # contains marginal cost for nonbasic columns in the tableau and there corresponding columns
	track = zeros((2,n-1)) # tracks nonbasic column in tableau.
			       # first row contains actual no of column and 2nd row contains it's position in table

	j = 0	
	for i in range(n):
		if(not(td[i] in bd)):
			
			track[0,j] = td[i]
			track[1,j] = i
			mcost[j] = dot(table[:,i],gb) -g[td[i]-1]
			if(td[i] in eqc):
				mcost[j+n-1] = -1*mcost[j]
			else:
				mcost[j+n-1] = 2*z-1*mcost[j]
			j = j+1	
	return mcost,track,z


def gj(table,i,j,m,n):
	"""gauss-jordan elimination"""
	
	table[i,:] = table[i,:]/table[i,j]
	for k in range(m+1):
		if (k!=i):
			table[k,:] = table[k,:] - table[i,:]*table[k,j]
			
def getresult(bd,track,mcost,z,m,n):
	"""calculating final result"""
	
	a = zeros(m+1)
	for i in range(m):
		if(i+n+1 in bd):
			a[i] = -z
		elif(i+n+1+m+n in bd):
			a[i] = z
		else:
			if (i+n+1 in track[0]):
				for k in range(n-1):
					if(track[0,k]==i+n+1):
						j = k
			elif (i+n+1+m+n in track[0]):
				for k in range(n-1):
					if (track[0,k]==i+n+1+m+n):
						j = k+n-1	
			a[i] = mcost[j]-z
	a[-1] = z	
	return a

def test1():
	"""
	Test function for solver.
	"""
	Aeq = zeros((2,4))
	beq = zeros(2)
	Aeq[0] = [7,-4,5,3]
	Aeq[1] = [-2,1,5,4]
	beq = array([-30,15])
	n,m = Aeq.shape
	w = solver(Aeq,beq)
	print w
    #print mat(Aeq)*(mat(w[:m]).transpose())

def test2():
    """
	Test function for solver.
	"""
    Aeq = zeros((10,11))
    beq = zeros(10)
    Aeq[0] = [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]
    Aeq[1] = [0.0,-1.41421356237,-1.41421356237,0.0,-1.41421356237,0.0,1.41421356237,0.0,0.0,-2.0,2.0]
    Aeq[2] = [0.0,1.41421356237,-1.41421356237,0.0,-1.41421356237,-2.82842712475,1.41421356237,2.0,0.0,2.0,0.0]
    Aeq[3] = [0.0,2.0,2.0,2.82842712475,2.0,-1.73185491414e-16,2.0,2.0,2.82842712475,-1.73185491414e-16,-2.0]
    Aeq[4] = [0.0,-2.0,2.0,0.0,2.0,-0.0,2.0,0.0,0.0,-4.0,0.0]
    Aeq[5] = [0.0,2.82842712475,-2.82842712475,0.0,-2.82842712475,4.89842541529e-16,2.82842712475,4.0,0.0,-3.46370982829e-16,-0.0]
    Aeq[6] = [0.0,-2.82842712475,-2.82842712475,0.0,-2.82842712475,-0.0,2.82842712475,0.0,0.0,3.46370982829e-16,-4.0]
    Aeq[7] = [0.0,2.0,2.0,0.0,2.0,0.0,2.0,0.0,0.0,4.0,4.0]
    Aeq[8] = [0.0,2.0,2.0,0.0,2.0,8.0,2.0,4.0,0.0,4.0,0.0]
    Aeq[9] = [0.0,4.0,4.0,8.0,4.0,2.99932144364e-32,4.0,4.0,8.0,2.99932144364e-32,4.0]
    
    beq = array([-4.5,1.41421356237,-0.585786437627,-6.82842712475,-1.33226762955e-15,-2.0,4.82842712475,-6.0,-10.0,-18.0])
    n,m = Aeq.shape
    print Aeq
    print beq
    w = solver(Aeq,beq)
    print w
    print mat(Aeq)*(mat(w[:m]).transpose())

if __name__=="__main__":
    test2()
